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Abstract. We present a computational approach able to extract the vibron-vibron coupling 
strength in clusters. Our approach is based on ab initio Bom-Oppenheimer molecular 
dynamics — a method that is inherently anharmonic — and a projection formalism able to 
deliver the individual vibron occupation numbers. From the Fourier transform of the vibron 
energy autocorrelation function we obtain the coupling strength of each vibron to the most 
strongly coupled vibronic states. We apply the method to a SiioHig cluster and obtain a blue 
shift of the Si-Si vibrational modes with transverse acoustic character and a red shift of the 
other vibrational modes with increasing temperature. We link this behavior to the bond length 
expansion and the varying sign of the Griineisen parameter We find vibron-vibron coupling 
strength up to 2.5 THz with a moderate increase of about 5% when increasing the temperature 
from 50 to 150 K. 
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1. Introduction 

Colloidal semiconductor nanoclusters (NCs) have undergone an enormous development in 
the fields of optoelectronics, spintronics, and bio-labeling over the past two decades. The 
good control of the NC size makes it possible to tailor their electronic and optical properties. 
Many successful applications in these fields were already reported [HI IH |3l IH |5l |6l |3. The 
experiments are usually performed at room temperature, making the treatment of vibrational 
properties, especial the anharmonic effects, crucial [|8] [TOl. A solid understanding 
of the anharmonic effects of vibrations, is therefore a decisive step for the real world 
application of nanostructures, where the physical properties such as thermal conductivity 
in nanowires [fTTI . non-radiative transition via phonon in NCs [[T2l . and Raman spectra 
broadening in nanostructures [iTSl are dominated by the phonon lifetime. There are basically 
two approaches to address this problem theoretically. 

The first one is to calculate the electron-phonon coupling [[T4l [TSl [T6l [TTl and 
phonon-phonon coupling [[TSl [T9l [20l |2T1 l22l| terms via perturbative approaches. The 
temperature effects are subsequently included using the Bose-Einstein distribution of the 
lattice vibrations (phonons). The physical properties, such as spectral broadening, spin- 
flip, loss of quantum coherence, relaxation of charge carriers, and thermal conductivity 
have been studied theoretically [[T9l l20l [T6l [TTl |24l, mostly considering harmonic or quasi- 
harmonic (harmonic approximation performed at different volumes) approximations. Some 
have considered the third-order and fourth-order derivatives of the total energy [|2T1 l22l . but 
the large computational demand for these higher-order derivatives limit the study to very small 
nanostructures [|25l or bulk materials [IITI l22l . The 2n + 1 theorem, on the other hand limits 
the studies based on density functional perturbation theory (DFPT) to third order process [|26l . 

The second one is to use molecular dynamics simulations ll27l l28l |29l [30l . where the 
temperature effects including all anharmonic effects are intrinsically contained in the atomic 
trajectories of the simulations. The temperature-dependent vibrational density of states (DOS) 
and the thermal conductivity of nanostructures have already been studied using classical 
molecular dynamics simulations [|48l . The required empirical force field potentials are 
limited by the lack of transferability to different systems and by the inability to correctly 
predict chemical bonding processes. These difficulties are overcome in the case of ab 
initio molecular dynamics (AIMD) simulations (281 HH [30l [SB EH |20l ESI, which does not 
build upon the harmonic approximation but allows for a self-consistent rearrangement of the 
changes, including all anharmonic effects. Based on the accurately calculated forces from the 
electronic structure calculations, AIMD enables us to monitor the changes in the electronic 
and vibrational properties on-the-fly, with thermal effects included directly. Although AIMD 
simulations have been successfully applied to study a variety of problems Il34l . and the phonon 
lifetime in bulk has been studied using classical molecular dynamics simulations Il35l[36ll48l . 
a vibron-vibron interaction extracted from AIMD simulation is still lacking. 

In this work, we perform AIMD simulations within the Bom-Oppenheimer (BO) 
approximation of a SiioHie cluster to study the geometry and the vibrational spectra from 
Fourier transformed velocity auto-correlation functions. The converged vibrational modes 
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Figure 1. SiioHig cluster (sila-adamantane). Silicon atoms are drawn as large yellow spheres 
and hydrogen atoms as small white spheres. 

are obtained from a trajectory of about 46 ps (corresponding to 96,000 steps). We find a 
blue-shift of the Si-Si vibrational modes with transverse acoustic (TA) characters and a red- 
shift of the other vibrational modes with increasing temperature. We also see a broadening 
of all the vibrational modes with increasing temperature. The former can be linked to the 
negative (blue-shift) and positive (red-shift) Griineisen parameters along with the extended 
bond lengths. The latter is attributed to the low symmetry (proximity to the surface) enhancing 
anharmonic effects at high temperatures. We further present a computational approach that 
enables the extraction of inter- vibron coupling strength, taking all the anharmonic effects into 
account. We find, for our cluster, couplings in the range of 0. 15 to 2.5 THz which correspond 
to "Rabi" oscillations of the vibron occupation number in the range of 0.4 to 6 ps. 

2. Theoretical method 

2.1. AIMD simulations 

We construct a sila-adamantane (SiioHie) cluster that has the point group symmetry [|38l 
[39ll . as shown in Figure \T\ In this cage-shaped cluster, four silicon atoms are bonded to 
three other silicon atoms and terminated by one hydrogen atom, while the remaining six 
silicon atoms are bonded to two silicon atoms and saturated with two hydrogen atoms. The 
simulations are performed using density functional theory (DFT) within the local density 
approximation (LDA) and Trouiller-Martin norm-conserving pseudopotentials with an energy 
cutoff of 20 Ry ll40l . The supercell is simple cubic with an extent of 15 A in each direction. 
The initial geometry of the cluster is optimized until the minimum forces acting on the Si and 
H atoms are less than 3x 10~^ a.u. under constrained symmetry. 

The AIMD simulations are initially performed at a constant temperature (NVT- 
ensemble) using the Nose-Hoover chain thermostat pT| with time step of 20 a.u. (about 
0.48 fs) in order to equilibrating the system. Following the 2 ps equilibration in the NVT- 
ensemble, we start a constant energy (NVE-ensemble) simulation and the trajectories are 
recorded at each time step. All the AIMD simulations are performed with the CPMD 
code [|40ll . In order to improve the statistics, we chop the NVE-ensemble simulation into 
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Figure 2. Schematic diagram of the AIMD simulation process. We start by an NVT 
simulation of 5.75 ps. From this simulation, we extract 15 initial configurations for the NVE 
runs. These are taken after an initial 2 ps of simulation time, with a time interval At of 0.25 
ps. We perform 15 NVE simulations ("slices") lasting 3 ps (12 ps) for the calculation of 
vibrational DOS (vibrational cooling). For each slice we perform a simple moving average, as 
sketched on the lower right. 
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several slices. Each slice starts from a different NVT-ensemble equilibration time and ends 
with the same number of NVE-ensemble simulation time steps. 

2.2. Vibrational DOS 

In order to obtain the temperature-dependent vibrational DOS from the AIMD simulation, we 
calculate the velocity auto-correlation function C{t) [j42l, 
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where, ( ) denotes the time averaged value calculated along the entire trajectory, Vi^t^j) is 
the velocity vector of atom i in slice k at time origin point j. The number of atoms in the 
cluster, the number of time origin points, and the number of slices are labeled as A^, rit, and 
Us, respectively. The time origin points labeled as Iqi, to2, toj and ton* are given in Figure |2] 
In the present work, the number of slices Us is taken to be 15 with 6400 time steps and 3200 
time origin points in each slice. The NVE-ensemble simulations are performed for a total of 
96000 time steps corresponding to about 46 ps simulation time. One slice corresponds to a 
simulation time of approximately 3 ps. 

The vibrational DOS are calculated using the Fourier transform [|42l 

1 /• + 00 
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Figure 3. Convergence of the vibrational modes AVDOSk{i^) with frequency between (a) 
450 and 650 cm~^ (bending Si-H and rotation H-Si-H modes), and (b) 1900 and 2100 cm"! 
(stretching Si-H modes) at T = 800 K as function of the simulation time given in unit of slices 
(3 ps). The vertical axis describes the deviation of the vibrational DOS between a simulation 
with k and a simulation with fc — 1 time slices. The red color (maximum deviations) reflects a 
deviation of 30 %. 



and for demonstration of the convergence with the number of slices we use the change of the 
vibrational DOS 

AVDOSkiu) = \VDOSk{uj) - VDOSk^i{uj)\ , (3) 

showing the deviation of the vibrational DOS between a simulation with k and a simulation 
with k — 1 time slices. Figure |3] (a) and (b) shows AVDOSk{uj) for the low (a) and high (b) 
frequency range at a temperature of 800 K. In this figure, the maximum deviation is 30%. We 
observe that the vibrational modes start to converge after 35 ps (12 slices). 

Since the limited statistic we obtained from AIMD at low temperature even for a 
small cluster, we use the DFPT results based on the harmonic approximation of lattice 
dynamics as a low temperature limit. In this case, the DFPT represents a very good 
classical (neglecting zero-point motion) approximation. The vibrational frequencies u and 
the vibrational eigenvectors Ui are obtained from the eigenvalue equation: P3l l44l 

where, M^q) is the mass of atom i (j), V is the potential, Ri(j) denotes the atomic 
positions. The dynamical matrix elements ^^^^^^ dR dR calculated using linear response 

as implemented in the CPMD code l|40ll . In this approach, the second derivative of potential 
is computed from the linear response of the system to an infinitesimal replacement based on 
the perturbation theory [|26l . 



2.3. Energy relaxation 

The extraction of the potential energy of a certain vibrational mode, Ep{t), was first proposed 
by Ladd et al. [[35l and subsequently modified by McGaughey et al. If36l to include the total 
energy of each mode instead of the potential energy only ||36l [371 . In the quasi-harmonic 
approximation, where the changes in bond length due to thermal expansion are included but 
further anharmonic effects are excluded [|20l , the total energy of each vibrational mode E''{t) 
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— proportional to the occupation number of the vibrational mode — is expressed in terms of 
the time-dependent normal coordinates Q„{t) 

E'it) = liQim^t) + ooiQimM, (5) 

where 

Q'^it) = T.yji^^-i^^t)-R^]. (6) 

i 

The three-component vectors Ul' in Eq. © represent the three components belonging to atom 
i of the vibrational eigenvectors obtained from Eq. (HJ. is the equilibrium position of atom 
i obtained from the trajectory of the AIMD simulation as, 

Based on the quasi-harmonic approximation, the vibrational vectors Ul' used in Eq. ^ are 
calculated using DFPT with the atomic positions obtained from Eq. ([7]), i.e. from an AIMD 
simulation at a certain temperature. 

The attenuation of the vibrational amplitude reflects the mode relaxation processes and 
can be described by the auto-correlation function of the energy fluctuation, written as 

where, ( ) denotes the time averaged value calculated along the entire trajectory and 
5E''{t) = E^{t) — E" is the deviation from the average vibrational energy E" . 



3. Temperature dependent vibrational properties 

We now describe the temperature dependence of the vibrational properties. In Figure IH we 
plot the converged vibrational DOS calculated using Eq. © and © at (a) 1600 K, (b) 1400 K, 
(c) 1200 K, (d) 1000 K, and (e) 800 K along with (f) the zero temperature linear-response 
results from Eq. dD (all ignore the zero point vibrations). In this work, we use the linear- 
response results as a low temperature limit because of the limited statistic we obtained from 
AIMD at low temperature even for the small cluster. The vibrational eigenmodes obtained 
from the linear-response calculations are analyzed in terms of their displacement pattern by 
a projection onto bulk phonon modes (see Ref. [|45l ) and are divided into acoustic -like Si-Si 
modes and optical-like Si-Si modes. The vibrations between the silicon and hydrogen atoms 
are divided into bending Si-H and rotation H-Si-H modes, shear H-Si-H modes and stretching 
Si-H modes according to the displacement of the atoms. These vibrational modes are listed in 
Table [Hand labeled as M1-M5 in Figure |4](f). In the AIMD simulations, the character of the 
vibrational eigenmodes cannot be obtained directly B6l . 

From Figure HI we observe that the vibrational DOS of all the modes, especially for the 
Si-H modes, show a broadening with increasing temperature and that the vibrational density of 
the shear H-Si-H modes and the stretching Si-H modes decreases with increasing temperature. 
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Table 1. Vibrational modes of SiipHig cluster. 





label 


frequency (cm ^) 


acoustic-like Si-Si modes 
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optical-like Si-Si modes 
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Figure 4. Vibrational DOS of a SiioHig cluster obtained from AIMD simulations at 
the temperatures of (a) 1600 K, (b) 1400 K, (c) 1200 K, (d) 1000 K, and (e) 800 K. The 
vibrational modes of a SiioHig cluster obtained from linear-response (LR) calculations at zero 
temperature is given in (f) including the mode assignment. 

There is also a red shift with increasing temperature of all the vibrations except for acoustic- 
like Si-Si modes, which correspond to the TA phonon modes of bulk Si. 

To explain the reason for this temperature dependence, we plot in Figure |5] the average 
bond lengths of (a) Si-Si bonds and (b) Si-H bonds as a function of the distance of the 
respective bond center to the cluster center at T = 800, 1200, and 1600 K. In the cage- 
shaped SiioHie cluster, the bond lengths of each Si-Si bond (Si-H bond) are uniform after 
geometry optimization at zero temperature. We plot these optimized bond lengths as dashed 
lines in Figure |5l We see from Figure |5] (a) that the Si-Si bond lengths increase with 
increasing temperature. Moreover, the cluster shows a large bond length distribution at high 



Vibron-vibron coupling from AIMD 



8 



2.45 



CD 

c 
_aj 

T3 

c 
o 

03 

> 
< 



2.40- 



2.35 



(a) Si-Si bond 
o 800 K 
□ 1200 K 
A 1600 K ^ 

A 



8 ^ 



So 



Q. O □ 

„ o 



(b) Si-H bond 

A 

AD A 

^□^ A°°8 

OOa?° 

-5-D----fe-a6-«a 



1.65 — 

1.60 ^ 
"a 
c 
o 

J2 
03 
O) 

1.50 I 



1.55 



2.1 2.2 2.3 2.4 3.1 3.2 3.3 

Distance to cluster center (A) 



3.4 



Figure 5. Thermally averaged bond length distribution of (a) Si-Si bond and (b) Si-H bond as 
a function of their distance (bond center) to the dot center at T = 800 K (black circles), 1200 K 
(red squares), and 1600 K (green triangles). The optimized bond lengths of the SiioHig cluster 
at T = K are given as dashed lines. 



temperatures. The increase of bond length and its large distribution along with the positive 
Griineisen parameters (describing the change in phonon frequencies with volume) for the 
optical branches and longitudinal acoustic (LA) branch result in the softening (red shifting) 
and broadening of the optical-like and the LA-like Si-Si modes with increasing temperature. 
In contrast to the relatively small extension of Si-Si bonds 3%), the Si-H bonds show a 
large extension and distribution with increasing temperature. A bond length extension of as 
much as 9% at T = 1600 K is seen in Figure [5] (b). We attribute the large extension of the 
Si-H bonds to the geometry of the SiioHig cluster. In contrast to the silicon atoms, which 
are localized at the center of the tetrahedral structure (see Figure [B, the hydrogen atoms at 
the cluster surface are free to move outwards. This introduces a potential asymmetry towards 
the vacuum side and increases the anharmonicity. A relatively large extension of Si-H bonds 
is therefore obtained at high temperature, which explains the red shift and the broadening of 
the Si-H vibrational modes, considering the positive Griineisen parameters. Especially, the 
vibrations of the high frequency H-Si-H shear and Si-H stretching modes mainly consist of 
the displacements of hydrogen atoms. This leads to a significant reduction and broadening of 
the vibrational peaks corresponding to the H-Si-H shear modes and Si-H stretching modes in 
Figure |4](a)-(e). Finally, the blue shift of the TA-like Si-Si modes observed in Figure |4]is the 
result of the bond length extension along with the negative Griineisen parameters of the TA 
modes. 



4. Vibron-Vibron Coupling 

The vibrational energy of a certain mode is proportional to the mode occupation number 
E'^{t) oc n'^{t) (see Eq.[5l), and thus the time evolution of E'^{t) carries the information about 
the coupling between different vibrational modes. In the following we will refer to E'^{t) 
as to the "occupation auto-correlation function", to avoid confusions with the many variants 



Vibron-vibron coupling from AIMD 



9 




c 5 10 5 10 




5 10 0.25 0.5 0.75 

Time (ps) 



Figure 6. Occupation autocorrelation functions of vibrational modes with frequencies (a) 
61.87 cm-\ (b) 158.85 cm-\ (c) 339.52 cm-\ (d) 516.23 cm-\ (e) 673.16 cm-\ and (f) 
2102.95 cm-i at T = 100 K. 

of the energy auto-correlation function. The total vibrational enegy, J2u ^ conserved in 
our NVE simulation and energy is allowed to flow between vibrational modes. To extract 
a meaningful statistical average out of E'^it) we use the correlation function C^(t) (Eq. [8]), 
which decays to zero for large times t, when the occupation number is unrelated to the initial 
situation (at time to)- We have calculated the occupation auto-correlation functions for all the 
vibrational modes of the SiioHie cluster at T = 50, 100, 125, and 150 K. In figure [6] (a)-(f) 
we plot the vibrational occupation autocorrelation function of six selected vibrational modes 
of the SiioHie cluster at T = 100 K, which will be discussed in more detail. 

It should be pointed out that the AIMD simulations in this work are performed only at 
relatively low temperatures in order not to violate the quasi-harmonic approximation which 
is used to estimate the vibrational eigenmode energy [|35ll36ll . Indeed, the vibrational vectors 
are calculated using DFPT with the average finite-temperature atomic positions obtained 
from Eq. (|7]). Imaginary vibrational frequencies are introduced when the atoms deviate far 
from their equilibrium positions at high temperatures. 

First, we have to prominently assert that a lifetime, as can be extracted from a damped 
Rabi oscillation or from transitions treated at the level of Fermi's golden rule does not exist 
in our cluster. The NVE simulation does not allow for energy dissipation, that would possibly 
lead to exponential decay of high energy phonon modes, from which a lifetime could be 
extracted. On the other hand, the discrete and energetically sparse nature of the vibrons does 
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not allow for a treatment following Fermi's golden rule, where transitions into a perfectly 
dissipative continuum are assumed. Our occupation numbers describe the time evolution of 
a many vibronic level system without dissipative term. The change of occupation of the 
individual levels reflects the interlevel coupling. The strongest coupling will dominate the 
short time evolution. 

In figure [6] (a)-(f) we observe a rather different behavior of the different vibrational 
modes. In all cases, however, there is an ultrafast oscillation superimposed on a slower 
variation. To analyze these results in a quantitative way, we replot in figure |7] (a) the 
occupation autocorrelation function for the vibrational mode with frequency 387.9 cm^^ 
(Fig. [6] (c)) along with its Fourier transform (panel b). The Fourier analysis shows that the 
high frequency oscillation (Period II in|7](a)), in this case at 11.86 THz, corresponds to the 
vibron frequency. The origin of this oscillation resides in errors in the projection, Eq. (|6j, has 
no physical meaning, and will not be discussed further. The slower variation (Period I in|7](a)) 
is very clear from Fig. |7](b) and has a frequency of 1.89 THz. We identify this feature as the 
frequency with which the vibron state undergoes periodic oscillations with strongly coupled 
vibronic states. It describes how the specific vibron mode can transfer energy to other modes, 
i.e., decay into lower energy vibrations or how two, or more, low energy vibrations can create 
a high energy vibron. Both processes are present in our AIMD calculations. The fact that we 
are dealing with a many-level (we have 78 vibrational eigenmodes) system explain that we 
do not see clear Rabi oscillations between two distinct levels but a fast oscillation (1.89 THz) 
modulated in time by the interaction with more weakly coupled levels. We cannot extract 
the weak coupling components from our results (that are all hidden in the low frequency 
Fourier transform Figure E](b), and would require much longer simulation times to be properly 
captured), but can clearly obtain the strongest coupling. 

In Figure |7] (c) we summarize our Fourier analysis of all the vibron modes and plot the 
oscillation period as a function of the vibron frequency. We find strongly coupled (short 
oscillation periods) modes at all available frequencies, starting from the bulk acoustic-like 
modes to the passivant hydrogen vibrations. The coupling strength varies between 0.15 and 
2.5 THz. 

We now study the temperature effects on the coupling strength of a certain vibrational 
mode. In Figure [8] (a)-(d), we plot the occupation autocorrelation functions of a certain 
acoustic-like mode with varying temperatures and in Fig. |8le) the extracted coupling 
strengths. We find that the vibrational frequencies show a red shift with increasing 
temperature which can be related to the thermal expansion and behaves as expected, and that 
the coupling strength increases slightly, from 2.06 to 2. 16 THz. This increase can be attributed 
to higher-order anharmonic effects [|47l l48l . With the increase in temperature, the effects of 
higher order anharmonicity, which includes four-vibron, five- vibron, and even higher order 
interactions, become stronger [|32l . 
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Figure 7. (a) Occupation autocorrelation function of one vibrational mode (see Figure |6] 
(c)) with labeled vibron-vibron oscillation period (Period I) and vibrational oscillation period 
(Period II). (b) Fourier transform of the occupation autocorrelation function of (a) vibrational 
mode with labeled Peak I corresponding to the vibron-vibron coupling strength and peak II 
coresponding to the vibrational frequency, (c) The vibron-vibron oscillation period for all the 
vibrational modes (circles) at 100 K along with the bulk phonon lifetime (dashed line) from 
experiment. 



5. Summary 

We performed first-principles AIMD simulations, within the BO approximation, to study 
the temperature dependent vibrational properties of a SiioHig cluster. We first calculate 
the vibrational DOS using the Fourier transformed velocity autocorrelation functions and 
compare the results to DFPT, obtaining good agreement. We quantify the softening and 
broadening, with increasing temperature, of the Si-Si LA-like modes, the Si-Si optical modes 
and the Si-H modes (especially for the H-Si-H shear modes and Si-H stretching modes) and 
analyze the results in terms of bond length variations. Subsequently, we suggest a method 
to calculate the vibron-vibron coupling strength, including all the anharmonic effects, using 
a Fourier transformation of the occupation autocorrelation function. For our Si cluster, we 
find a coupling strength between 0.15 and 2.5 THz, which corresponds to oscillations in 
the occupation of vibron states (akin Rabi oscillations) in the range of 0.4 to 6 ps. The 
obtained coupling parameters could be used in further study of the dynamical processes in 
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Figure 8. Occupation autocorrelation functions of a certain vibrational mode (see Figure. |6] 
(b)) at temperature (a) 50 K, (b) 100 K, (c) 125 K, and (d) 150K. The vibron-vibron coupling 
strength obtained from the occupation autocorrelation functions are given in (e). 



nanosystem based on rate equations and potentially including dissipation through the coupling 
to an external phonon bath [|49l |50l. We conclude that, although our approach enables the 
extraction of anharmonic vibron coupling parameters otherwise not available, it still requires 
trajectories of around 24 ps, which presently limits it's applicability to rather small structures. 
Further challenges involve the non-adiabatic coupling to electronic states and the coupling to 
a dissipative environment, which we have ignored. 
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